import pandas as pd
import numpy as np
# import matplotlib.pyplot as plt
from util_set_zh_matplot import plt
import seaborn as sns
from pathlib import Path
import pdb
import os

import util_for_output_zh

# 在查看数据前设置
pd.set_option('display.max_columns', None)
pd.set_option('display.max_rows', None)
pd.set_option('display.width', None)
pd.set_option('display.max_colwidth', None)

def main():
    # stage1()
    # stage2()
    # stage3()
    # stage4()
    # stage5()
    # stage6()
    # stage7()
    # stage8()
    # stage8_1()
    # stage9()
    # stage9_1()
    # stage9_2()
    stage_9_2_optim()


def stage1():
    '''
    1, 读取 附件.xlsx , 读取其中的 两个sheet 保存csv 到 ques0_prepare_data 文件夹
    '''
    # 读取 附件.xlsx 文件
    excel_path = '附件.xlsx'
    output_dir = Path('ques0_prepare_data')
    output_dir.mkdir(exist_ok=True)

    # 获取所有sheet名称
    xls = pd.ExcelFile(excel_path)
    sheet_names = xls.sheet_names

    # 只读取前两个sheet
    for i, sheet_name in enumerate(sheet_names[:2]):
        df = pd.read_excel(excel_path, sheet_name=sheet_name)
        csv_path = output_dir / f'sheet{i+1}_{sheet_name}.csv'
        df.to_csv(csv_path, index=False, encoding='utf-8-sig')

'''
两张表
ques0_prepare_data/sheet1_男胎检测数据.csv
ques0_prepare_data/sheet2_女胎检测数据.csv

head info
   序号  孕妇代码  年龄     身高    体重                 末次月经       IVF妊娠   检测日期  检测抽血次数   检测孕周      孕妇BMI    原始读段数  在参考基因组上比对的比例   重复读段的比例  唯一比对的读段数        GC含量  13号染色体的Z值  18号染色体的Z值  21号染色体的Z值   X染色体的Z值   Y染色体的Z值    Y染色体浓度    X染色体浓度  13号染色体的GC含量  18号染色体的GC含量  21号染色体的GC含量  被过滤掉读段数的比例 染色体的非整倍体 怀孕次数  生产次数 胎儿是否健康
0   1     A001    31      160.0   72.0          2023-02-01 00:00:00  自然受孕  20230429       1  11w+6  28.125000  5040534      0.806726  0.027603     3845411  0.399262   0.782097  -2.321212  -1.026003 -0.062103 -1.035610  0.025936  0.038061     0.377069     0.389803     0.399399    0.027484      NaN    1     0      是
1   2  A001  31  160.0  73.0  2023-02-01 00:00:00  自然受孕  20230531       2  15w+6  28.515625  3198810      0.806393  0.028271     2457402  0.393299   0.692856   1.168521  -2.595099  0.582183 -0.363519  0.034887  0.059572     0.371542     0.384771     0.391706    0.019617      NaN    1     0      是

'''
# 处理男胎数据
def stage2():
    # 处理男胎数据
    # 1 输出 存在 "缺失值" 列，每个列 存在多少个 缺失值
    # 读取男胎数据
    male_data_path = Path('ques0_prepare_data/sheet1_男胎检测数据.csv')
    df_male = pd.read_csv(male_data_path, encoding='utf-8-sig')

    # 统计每一列的缺失值数量
    missing_counts = df_male.isnull().sum()

    # 只输出存在缺失值的列及其缺失数量
    missing_counts = missing_counts[missing_counts > 0]
    if not missing_counts.empty:
        print("存在缺失值的列及缺失数量：")
        print(missing_counts)
    else:
        print("男胎数据中没有缺失值。")
    
'''
存在缺失值的列及缺失数量：
末次月经         12
染色体的非整倍体    956
dtype: int64
'''
def stage3():
    # 末次月经 缺失值NAN 补全为 末次月经日期 = 检测日期 - 孕周天数
    #  df_filled.loc[missing_lmp, '末次月经'] = df_filled.loc[missing_lmp].apply(
        # lambda row: calculate_last_menstrual_date(row['检测日期'], row['孕周天数']), axis=1)
    # df['染色体的非整倍体'] = df['染色体的非整倍体'].fillna("正常")   
    # 补全缺失值逻辑
    from util_for_ques0 import parse_gestational_age, calculate_last_menstrual_date

    male_data_path = Path('ques0_prepare_data/sheet1_男胎检测数据.csv')
    filename_fillnan = 'sheet1_男胎检测数据_fillnan'
    df_male = pd.read_csv(male_data_path, encoding='utf-8-sig')

    # 1. 解析孕周为天数
    if '检测孕周' in df_male.columns:
        df_male['孕周天数'] = df_male['检测孕周'].apply(parse_gestational_age)
    else:
        df_male['孕周天数'] = None

    # 2. 末次月经缺失值补全
    if '末次月经' in df_male.columns and '检测日期' in df_male.columns and '孕周天数' in df_male.columns:
        missing_lmp = df_male['末次月经'].isnull()
        if missing_lmp.sum() > 0:
            df_male.loc[missing_lmp, '末次月经'] = df_male.loc[missing_lmp].apply(
                lambda row: calculate_last_menstrual_date(row['检测日期'], row['孕周天数']), axis=1
            )
            # 若仍有缺失，使用众数填充
            if df_male['末次月经'].isnull().sum() > 0:
                lmp_mode = df_male['末次月经'].mode()[0]
                df_male['末次月经'].fillna(lmp_mode, inplace=True)

    # 3. 染色体的非整倍体缺失值补全为"正常"
    if '染色体的非整倍体' in df_male.columns:
        df_male['染色体的非整倍体'] = df_male['染色体的非整倍体'].fillna("正常")

    # 4. 删除临时列
    if '孕周天数' in df_male.columns:
        df_male.drop('孕周天数', axis=1, inplace=True)

    # 5. 保存补全后的数据
    df_male.to_csv(f'ques0_prepare_data/{filename_fillnan}.csv', index=False, encoding='utf-8-sig')
    print(f"补全后的男胎数据已保存到 ques0_prepare_data/{filename_fillnan}.csv")

    # 检查补全后是否还存在缺失值，并输出
    missing_info = df_male.isnull().sum()
    missing_cols = missing_info[missing_info > 0]
    if not missing_cols.empty:
        print("补全后仍存在缺失值的列及数量：")
        print(missing_cols)
    else:
        print("补全后数据中已不存在缺失值。")

'''
补全后的男胎数据已保存到 ques0_prepare_data/sheet1_男胎检测数据_fillnan.csv
补全后数据中已不存在缺失值。
'''

# 处理女胎数据
def stage4():
    # 处理女胎数据
    # 1 输出 存在 "缺失值" 列，每个列 存在多少个 缺失值
    # 假设女胎数据文件路径
    female_data_path = 'ques0_prepare_data/sheet2_女胎检测数据.csv'
    try:
        df_female = pd.read_csv(female_data_path, encoding='utf-8-sig')
    except Exception as e:
        print(f"读取女胎数据文件失败: {e}")
        return

    # 统计每一列的缺失值数量
    missing_info = df_female.isnull().sum()
    missing_cols = missing_info[missing_info > 0]
    if not missing_cols.empty:
        print("存在缺失值的列及其缺失数量：")
        print(missing_cols)
    else:
        print("女胎数据中不存在缺失值。")

'''
存在缺失值的列及其缺失数量：
孕妇BMI            1
Unnamed: 20    605
Unnamed: 21    605
染色体的非整倍体       538
dtype: int64
'''
def stage5():
    # 处理缺失值
    # 删除Unnamed列 BMI 缺失较少，删除不影响整体，所以删除缺失的BMI所在行
    # 染色体的非整倍体 像 上面一样补全 为 正常 df['染色体的非整倍体'] = df['染色体的非整倍体'].fillna("正常") 
    # 读取女胎数据
    female_data_path = 'ques0_prepare_data/sheet2_女胎检测数据.csv'
    try:
        df_female = pd.read_csv(female_data_path, encoding='utf-8-sig')
    except Exception as e:
        print(f"读取女胎数据文件失败: {e}")
        return

    # 删除所有列名中包含 'Unnamed' 的列
    unnamed_cols = [col for col in df_female.columns if 'Unnamed' in col]
    if unnamed_cols:
        print(f"删除列: {unnamed_cols}")
        df_female.drop(columns=unnamed_cols, inplace=True)

    # 删除孕妇BMI缺失的行
    if '孕妇BMI' in df_female.columns:
        bmi_missing = df_female['孕妇BMI'].isnull().sum()
        if bmi_missing > 0:
            print(f"删除孕妇BMI缺失的 {bmi_missing} 行")
            df_female = df_female[~df_female['孕妇BMI'].isnull()]

    # 补全“染色体的非整倍体”缺失值为“正常”
    if '染色体的非整倍体' in df_female.columns:
        missing_chrom = df_female['染色体的非整倍体'].isnull().sum()
        if missing_chrom > 0:
            print(f"补全染色体的非整倍体缺失值 {missing_chrom} 个为 '正常'")
            df_female['染色体的非整倍体'] = df_female['染色体的非整倍体'].fillna("正常")

    # 检查补全后是否还存在缺失值，并输出
    missing_info = df_female.isnull().sum()
    missing_cols = missing_info[missing_info > 0]
    if not missing_cols.empty:
        print("补全后仍存在缺失值的列及数量：")
        print(missing_cols)
    else:
        print("补全后数据中已不存在缺失值。")

    # 保存处理后的女胎数据
    output_path = 'ques0_prepare_data/sheet2_女胎检测数据_fillnan.csv'
    df_female.to_csv(output_path, index=False, encoding='utf-8-sig')
    print(f"补全后的女胎数据已保存到 {output_path}")

'''
删除列: ['Unnamed: 20', 'Unnamed: 21']
删除孕妇BMI缺失的 1 行
补全染色体的非整倍体缺失值 537 个为 '正常'
补全后数据中已不存在缺失值。
补全后的女胎数据已保存到 ques0_prepare_data/sheet2_女胎检测数据_fillnan.csv
'''

def stage6():
    # 对 补全后的两张表
    # 按下述规则 初步筛选，筛选后的表 存储为 {fill_name}_sifted.csv
    '''
    原始读段数 < 2000000
    孕周数 10~25w 之间
    '''
    # 读取补全后的男胎和女胎数据
    male_path = 'ques0_prepare_data/sheet1_男胎检测数据_fillnan.csv'
    female_path = 'ques0_prepare_data/sheet2_女胎检测数据_fillnan.csv'
    df_male = pd.read_csv(male_path, encoding='utf-8-sig')
    df_female = pd.read_csv(female_path, encoding='utf-8-sig')

    # 定义筛选条件
    # 原始读段数 < 2000000
    # 孕周数 10~25w 之间
    def parse_weeks(weeks_str):
        # 解析孕周字符串如 '11w+6' 为周数（浮点数）
        import re
        if pd.isna(weeks_str):
            return None
        match = re.match(r'(\d+)w(?:\+(\d+))?', str(weeks_str))
        if match:
            weeks = int(match.group(1))
            days = int(match.group(2)) if match.group(2) else 0
            return weeks + days / 7
        return None

    # 对男胎和女胎分别处理
    for df, name in zip([df_male, df_female], ['男胎', '女胎']):
        # 解析孕周为周数
        df['孕周_周'] = df['检测孕周'].apply(parse_weeks)
        before_shape = df.shape

        # 筛选
        df_sifted = df[
            (df['原始读段数'] >= 2000000) &
            (df['孕周_周'] >= 10) &
            (df['孕周_周'] <= 25)
        ].copy()

        after_shape = df_sifted.shape
        print(f"{name}数据: 原始{before_shape}，筛选后{after_shape}")

        # 删除临时列
        df_sifted.drop(columns=['孕周_周'], inplace=True)

        # 保存
        out_path = f'ques0_prepare_data/sheet{"1" if name=="男胎" else "2"}_{name}检测数据_fillnan_sifted.csv'
        df_sifted.to_csv(out_path, index=False, encoding='utf-8-sig')
        print(f"{name}筛选后数据已保存到 {out_path}")

'''
男胎数据: 原始(1082, 32)，筛选后(1067, 32)
男胎筛选后数据已保存到 ques0_prepare_data/sheet1_男胎检测数据_fillnan_sifted.csv
女胎数据: 原始(604, 30)，筛选后(593, 30)
女胎筛选后数据已保存到 ques0_prepare_data/sheet2_女胎检测数据_fillnan_sifted.csv
'''

def stage7():
    # 处理 男胎数据 在 Y染色体浓度 和 孕妇BMI 上的 异常值
    # 根据stage3_outlier_process函数 中的 z-score 方法处理异常值, 处理好的数据 存放在 ques1_prepare_data/{filename}_outlier.csv
    # Y染色体浓度 箱线图 和 孕妇BMI 箱线图 保存到 ques1_prepare_data 文件夹下
    # X轴为孕周 Y轴为 Y染色体浓度 and 孕妇BMI 绘制 双折线图 ques1_prepare_data 文件夹下
    import os

    def _stage7():
        """
        处理男胎数据在 Y染色体浓度 和 孕妇BMI 上的异常值
        - 使用z-score方法处理异常值
        - 处理好的数据存放在 ques1_prepare_data/{filename}_outlier.csv
        - Y染色体浓度箱线图和孕妇BMI箱线图保存到 ques1_prepare_data 文件夹下
        - X轴为孕周，Y轴为 Y染色体浓度 和 孕妇BMI 绘制双折线图，保存到 ques1_prepare_data 文件夹下
        """
        import pandas as pd
        import numpy as np
        import matplotlib.pyplot as plt

        # 创建输出文件夹
        out_dir = "ques1_prepare_data"
        os.makedirs(out_dir, exist_ok=True)

        # 读取男胎数据
        infile = "ques0_prepare_data/sheet1_男胎检测数据_fillnan_sifted.csv"
        df = pd.read_csv(infile, encoding='utf-8-sig')

        # 只保留需要的列
        cols_needed = ['孕妇BMI', 'Y染色体浓度', '检测孕周']
        for col in cols_needed:
            if col not in df.columns:
                raise ValueError(f"缺少必要列: {col}")

        # 处理异常值（z-score法）
        def remove_outliers_zscore(df, col, z_thresh=3):
            vals = df[col]
            mean = vals.mean()
            std = vals.std()
            z_scores = (vals - mean) / std
            mask = z_scores.abs() <= z_thresh
            return df[mask].copy(), mask

        # 处理Y染色体浓度异常值
        df_y, mask_y = remove_outliers_zscore(df, 'Y染色体浓度', z_thresh=3)
        # 处理孕妇BMI异常值
        df_bmi, mask_bmi = remove_outliers_zscore(df, '孕妇BMI', z_thresh=3)

        # 取交集（两个都不是异常值的行）
        mask_both = mask_y & mask_bmi
        df_outlier = df[mask_both].copy()

        # 保存处理后的数据
        out_csv = os.path.join(out_dir, "sheet1_男胎检测数据_fillnan_sifted_outlier.csv")
        df_outlier.to_csv(out_csv, index=False, encoding='utf-8-sig')
        print(f"去除异常值后的男胎数据已保存到 {out_csv}")

        # 绘制箱线图
        plt.figure(figsize=(8, 6))
        plt.boxplot(df['Y染色体浓度'].dropna(), vert=True, patch_artist=True)
        plt.title("Y染色体浓度箱线图")
        plt.ylabel("Y染色体浓度")
        plt.savefig(os.path.join(out_dir, "Y染色体浓度_箱线图.png"), bbox_inches='tight')
        plt.close()

        plt.figure(figsize=(8, 6))
        plt.boxplot(df['孕妇BMI'].dropna(), vert=True, patch_artist=True)
        plt.title("孕妇BMI箱线图")
        plt.ylabel("孕妇BMI")
        plt.savefig(os.path.join(out_dir, "孕妇BMI_箱线图.png"), bbox_inches='tight')
        plt.close()

        # 解析孕周为周数
        def parse_weeks(weeks_str):
            import re
            if pd.isna(weeks_str):
                return np.nan
            match = re.match(r'(\d+)w(?:\+(\d+))?', str(weeks_str))
            if match:
                weeks = int(match.group(1))
                days = int(match.group(2)) if match.group(2) else 0
                return weeks + days / 7
            return np.nan

        df_outlier['孕周_周'] = df_outlier['检测孕周'].apply(parse_weeks)

        # 按孕周排序
        df_outlier = df_outlier.sort_values('孕周_周')

        # 绘制双折线图
        # 双折线图,左侧是Y染色体的量度,右边是BMI的量度
        # 并且显示平滑后的 两个折现
        # 重新绘制双折线图，左轴为Y染色体浓度，右轴为孕妇BMI，并对两条曲线做平滑处理
        # 如果同一孕周出现多个值,则取均值,其他值作为(min,max)呈现一个色域表示波动范围
        # 1. 按孕周分组，计算每组的均值、最小值、最大值
        grouped = df_outlier.groupby('孕周_周')
        x_unique = np.array(sorted(grouped.groups.keys()))
        y_y_mean = grouped['Y染色体浓度'].mean().values
        y_y_min = grouped['Y染色体浓度'].min().values
        y_y_max = grouped['Y染色体浓度'].max().values
        y_bmi_mean = grouped['孕妇BMI'].mean().values
        y_bmi_min = grouped['孕妇BMI'].min().values
        y_bmi_max = grouped['孕妇BMI'].max().values

        # 2. 对均值序列做平滑
        from scipy.ndimage import gaussian_filter1d
        y_y_mean_smooth = gaussian_filter1d(y_y_mean, sigma=1)
        y_bmi_mean_smooth = gaussian_filter1d(y_bmi_mean, sigma=1)

        # 3. 绘制双折线图，带色域
        fig, ax1 = plt.subplots(figsize=(20, 6))
        color1 = 'tab:blue'
        color2 = 'tab:red'

        # Y染色体浓度
        ax1.set_xlabel("孕周（周）")
        ax1.set_ylabel("Y染色体浓度", color=color1)
        # 原始点
        # ax1.scatter(df_outlier['孕周_周'], df_outlier['Y染色体浓度'], color=color1, alpha=0.3, label='Y染色体浓度原始')
        # 均值线
        l1 = ax1.plot(x_unique, y_y_mean, 'o-', color=color1, label='Y染色体浓度均值')
        # 平滑线
        l2 = ax1.plot(x_unique, y_y_mean_smooth, '-', color=color1, linewidth=2, label='Y染色体浓度平滑')
        # 色域
        ax1.fill_between(x_unique, y_y_min, y_y_max, color=color1, alpha=0.15, label='Y染色体浓度波动范围')
        ax1.tick_params(axis='y', labelcolor=color1)

        # 孕妇BMI
        ax2 = ax1.twinx()
        ax2.set_ylabel("孕妇BMI", color=color2)
        # 原始点
        # ax2.scatter(df_outlier['孕周_周'], df_outlier['孕妇BMI'], color=color2, alpha=0.3, marker='s', label='孕妇BMI原始')
        # 均值线
        l3 = ax2.plot(x_unique, y_bmi_mean, 's-', color=color2, label='孕妇BMI均值')
        # 平滑线
        l4 = ax2.plot(x_unique, y_bmi_mean_smooth, '-', color=color2, linewidth=2, label='孕妇BMI平滑')
        # 色域
        ax2.fill_between(x_unique, y_bmi_min, y_bmi_max, color=color2, alpha=0.15, label='孕妇BMI波动范围')
        ax2.tick_params(axis='y', labelcolor=color2)

        # 合并图例
        # 需要手动合并两个y轴的图例
        handles1, labels1 = ax1.get_legend_handles_labels()
        handles2, labels2 = ax2.get_legend_handles_labels()
        ax1.legend(handles1 + handles2, labels1 + labels2, loc='upper left')

        plt.title("孕周- Y染色体浓度 & 孕妇BMI 双折线图（双Y轴，均值+平滑+波动范围）")
        plt.tight_layout()
        plt.savefig(os.path.join(out_dir, "孕周_Y染色体浓度_孕妇BMI_双Y轴_平滑_双折线图.png"), bbox_inches='tight')
        plt.close()
    _stage7()

# 去除异常值后的男胎数据已保存到 ques1_prepare_data/sheet1_男胎检测数据_fillnan_sifted_outlier.csv
'''
题目如下
题2  临床证明，男胎孕妇的BMI是影响胎儿Y染色体浓度的最早达标时间（即浓度达到或超
过4%的最早时间）的主要因素。试对男胎孕妇的BMI进行合理分组，给出每组的BMI区间和最佳NIPT
时点，使得孕妇可能的潜在风险最小，并分析检测误差对结果的影响。
'''
# 使用题目假定的 BMI 分组 做的 数据处理
# 如何 自己根据数据 划分BMI 分组 -> 聚类
# 然后根据 聚类的 结果 就是 BMI 的 最终分组
# 聚类的中心是什么?
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
from scipy import stats
def stage8():
    # 创建输出文件夹
    out_dir = "ques2_prepare_data"
    os.makedirs(out_dir, exist_ok=True)
    # 先读取 sheet1_男胎检测数据_fillnan_sifted_outlier.csv 打印head info
    data = pd.read_csv('ques1_prepare_data/sheet1_男胎检测数据_fillnan_sifted_outlier.csv')
    # 检查缺失值情况
    missing_counts = data.isnull().sum()
    missing_cols = missing_counts[missing_counts > 0]
    if not missing_cols.empty:
        print("存在缺失值的列及缺失数量：")
        print(missing_cols)
    else:
        print("数据中没有缺失值。")

    '''
    时间统一
        末次月经 两种日期格式
        2023-02-01 00:00:00 -> 2023-02-01
        2023-02-25 -> 2023-02-25
        检测日期 两种日期格式
        2023-06-07 00:00:00 -> 2023-06-07
        20230702 -> 2023-07-02
    分别依次处理这里的日期
    '''
    # 处理'末次月经'列
    if '末次月经' in data.columns:
        def clean_lmp(x):
            if pd.isnull(x):
                return x
            x = str(x)
            # 2023-02-01 00:00:00 -> 2023-02-01
            if ' ' in x:
                x = x.split(' ')[0]
            return x
        data['末次月经'] = data['末次月经'].apply(clean_lmp)

    # 处理'检测日期'列
    if '检测日期' in data.columns:
        def clean_test_date(x):
            if pd.isnull(x):
                return x
            x = str(x)
            # 2023-06-07 00:00:00 -> 2023-06-07
            if ' ' in x:
                x = x.split(' ')[0]
            # 20230702 -> 2023-07-02
            if len(x) == 8 and x.isdigit():
                try:
                    x = f"{x[:4]}-{x[4:6]}-{x[6:8]}"
                except Exception:
                    pass
            return x
        data['检测日期'] = data['检测日期'].apply(clean_test_date)

    # pdb.set_trace()

    # 筛选男胎数据（Y染色体浓度>0.04即4%）
    data = data[data['Y染色体浓度'] >= 0.04].copy()

    # 检查日期列的数据类型，如果不是datetime64，则手动转换
    if data['末次月经'].dtype != 'datetime64[ns]':
        data['末次月经'] = pd.to_datetime(data['末次月经'], errors='coerce')  # errors='coerce'将无效格式转为NaT

    if data['检测日期'].dtype != 'datetime64[ns]':
        data['检测日期'] = pd.to_datetime(data['检测日期'], errors='coerce')
    
    # 清除日期无效的行
    # male_data = male_data.dropna(subset=['末次月经', '检测日期'])



    # 计算孕周（从末次月经到检测日期的周数）
    data['实际孕周'] = (data['检测日期'] - data['末次月经']).dt.days / 7
    male_data = data

    ## 2
    # 使用肘部法则确定最佳聚类数量
    bmi_values = male_data[['孕妇BMI']].values
    distortions = []
    K = range(1, 6)
    for k in K:
        kmeans = KMeans(n_clusters=k, random_state=42)
        kmeans.fit(bmi_values)
        distortions.append(kmeans.inertia_)
    
    # 绘制肘部法则图
    plt.figure(figsize=(10, 5))
    plt.plot(K, distortions, 'bx-')
    plt.xlabel('聚类数量')
    plt.ylabel('SSE (Sum of Squared Errors)')
    plt.title('肘部法则确定最佳聚类数量')
    plt.show()
    
    # 根据肘部法则选择k=3进行聚类
    kmeans = KMeans(n_clusters=3, random_state=42)
    male_data['BMI分组'] = kmeans.fit_predict(male_data[['孕妇BMI']])
    
    # 获取各组的BMI范围
    bmi_ranges = male_data.groupby('BMI分组')['孕妇BMI'].agg(['min', 'max'])
    bmi_ranges['区间'] = bmi_ranges.apply(lambda x: f"{x['min']:.1f}-{x['max']:.1f}", axis=1)
    
    # 重命名分组
    group_names = {
        0: "低BMI组",
        1: "中BMI组", 
        2: "高BMI组"
    }
    male_data['BMI分组'] = male_data['BMI分组'].map(group_names)

    print("\nBMI分组结果:")
    print(bmi_ranges[['区间']])
    plt.figure(figsize=(10, 6))
    sns.boxplot(data=male_data, x='BMI分组', y='孕妇BMI')
    plt.title('BMI分组分布')
    # plt.show()
    plt.savefig( os.path.join( out_dir , 'BMI分组分布.png' ) )
    plt.close()

    ## 3
    results = []
    
    for group in group_names.values():
        group_data = male_data[male_data['BMI分组'] == group]
        
        # 找到Y浓度首次≥4%的最小孕周
        first_reach = group_data.groupby('孕妇代码').apply(
            lambda x: x[x['Y染色体浓度'] >= 0.04]['实际孕周'].min()
        ).mean()
        
        # 计算95%置信区间
        ci_low, ci_high = stats.t.interval(0.95, len(group_data)-1,
                                          loc=first_reach,
                                          scale=stats.sem(group_data['实际孕周']))
        
        results.append({
            'BMI分组': group,
            '最佳NIPT孕周(均值)': round(first_reach, 1),
            '95%置信区间': f"{ci_low:.1f}-{ci_high:.1f}",
            '建议检测时间': f"{int(first_reach)}w+{int((first_reach%1)*7)}d"
        })

    # 4. 确定最佳NIPT时点
    time_df = pd.DataFrame(results)
    print("\n最佳NIPT时点建议:")
    print(time_df)

    # 5. 误差影响分析
    def error_analysis(male_data):
        # 模拟检测误差（±10%）
        np.random.seed(42)
        male_data['Y浓度误差'] = male_data['Y染色体浓度'] * np.random.uniform(0.9, 1.1, len(male_data))
        
        # 计算误差影响
        error_effect = male_data.groupby('BMI分组').apply(
            lambda x: (x['Y浓度误差'] >= 0.04).sum() / len(x)
        )
        
        return error_effect

    error_effect = error_analysis(male_data)
    print("\n检测误差影响（各组正确识别率）:")
    print(error_effect)

    # 6. 可视化Y浓度达标时间
    plt.figure(figsize=(12, 6))
    sns.lineplot(data=male_data, x='实际孕周', y='Y染色体浓度', 
                hue='BMI分组', ci=95)
    # 设置x轴范围为10到25
    plt.xlim(10, 25)
    plt.axhline(y=0.04, color='r', linestyle='--', label='4%阈值')
    plt.title('不同BMI组Y染色体浓度随孕周变化')
    plt.legend()
    # plt.show()
    plt.savefig( os.path.join( out_dir , '不同BMI组Y染色体浓度随孕周变化.png' ) )
    plt.close()

    
def stage8_1():
    data = pd.read_csv('ques1_prepare_data/sheet1_男胎检测数据_fillnan_sifted_outlier.csv')
    # 检查缺失值情况
    missing_counts = data.isnull().sum()
    missing_cols = missing_counts[missing_counts > 0]
    if not missing_cols.empty:
        print("存在缺失值的列及缺失数量：")
        print(missing_cols)
    else:
        print("数据中没有缺失值。")

    '''
    时间统一
        末次月经 两种日期格式
        2023-02-01 00:00:00 -> 2023-02-01
        2023-02-25 -> 2023-02-25
        检测日期 两种日期格式
        2023-06-07 00:00:00 -> 2023-06-07
        20230702 -> 2023-07-02
    分别依次处理这里的日期
    '''
    # 处理'末次月经'列
    if '末次月经' in data.columns:
        def clean_lmp(x):
            if pd.isnull(x):
                return x
            x = str(x)
            # 2023-02-01 00:00:00 -> 2023-02-01
            if ' ' in x:
                x = x.split(' ')[0]
            return x
        data['末次月经'] = data['末次月经'].apply(clean_lmp)

    # 处理'检测日期'列
    if '检测日期' in data.columns:
        def clean_test_date(x):
            if pd.isnull(x):
                return x
            x = str(x)
            # 2023-06-07 00:00:00 -> 2023-06-07
            if ' ' in x:
                x = x.split(' ')[0]
            # 20230702 -> 2023-07-02
            if len(x) == 8 and x.isdigit():
                try:
                    x = f"{x[:4]}-{x[4:6]}-{x[6:8]}"
                except Exception:
                    pass
            return x
        data['检测日期'] = data['检测日期'].apply(clean_test_date)

    print(data.head())
    # 统一日期后存回原来的地方
    data.to_csv('ques1_prepare_data/sheet1_男胎检测数据_fillnan_sifted_outlier.csv')

'''
题3  男胎Y染色体浓度达标时间受多种因素(身高、体重、年龄等)的影响，试综合考虑这些因
素、检测误差和胎儿的Y染色体浓度达标比例（即浓度达到或超过4%的比例），根据男胎孕妇的BMI，
给出合理分组以及每组的最佳NIPT时点，使得孕妇潜在风险最小，并分析检测误差对结果的影响。

'''
from sklearn.tree import DecisionTreeRegressor, plot_tree
def stage9():
    # 创建输出文件夹
    out_dir = "ques3_prepare_data"
    os.makedirs(out_dir, exist_ok=True)

    df = pd.read_csv('ques1_prepare_data/sheet1_男胎检测数据_fillnan_sifted_outlier.csv')
    # 1. 数据预处理
    # 假设df是已加载的DataFrame
    # 筛选男胎数据（Y染色体浓度>0）
    male_df = df[df['Y染色体浓度'] > 0].copy()

    # 计算孕周（从末次月经到检测日期）
    male_df['末次月经'] = pd.to_datetime(male_df['末次月经'])
    male_df['检测日期'] = pd.to_datetime(male_df['检测日期'])
    male_df['实际孕周'] = (male_df['检测日期'] - male_df['末次月经']).dt.days / 7

    # 2. 寻找最佳BMI分组
    # 先根据Y染色体浓度是否达标（>=0.04）生成"达标"列
    male_df['达标'] = (male_df['Y染色体浓度'] >= 0.04).astype(int)
    # 使用决策树寻找最佳分割点
    X = male_df[['孕妇BMI']].values
    y = male_df['达标'].values

    # 限制树深度为2，得到2-4个分组
    tree = DecisionTreeRegressor(max_depth=2, min_samples_leaf=10)
    tree.fit(X, y)

    # 提取分割点
    thresholds = np.sort(tree.tree_.threshold[tree.tree_.threshold != -2])
    print(f"自动计算的BMI分割点: {thresholds}")

    # 3. 分析每个BMI分组的达标情况
    bmi_bins = [0] + list(thresholds) + [np.inf]
    bmi_labels = [f"BMI≤{thresholds[0]:.1f}"] + \
                [f"{thresholds[i]:.1f}<BMI≤{thresholds[i+1]:.1f}" for i in range(len(thresholds)-1)] + \
                [f"BMI>{thresholds[-1]:.1f}"]

    male_df['BMI分组'] = pd.cut(male_df['孕妇BMI'], bins=bmi_bins, labels=bmi_labels)

    # 计算每个BMI分组在不同孕周的达标比例
    results = male_df.groupby(['BMI分组', '实际孕周'])['达标'].agg(['mean', 'count']).reset_index()
    results.rename(columns={'mean': '达标比例', 'count': '样本量'}, inplace=True)

    # 4. 确定每个BMI分组的最佳检测时间
    def find_optimal_week(group):
        # 优先选择达标比例≥80%的最早孕周
        candidates = group[group['达标比例'] >= 0.8]
        if not candidates.empty:
            return candidates['实际孕周'].min()
        
        # 如果没有达标比例≥80%的，选择达标比例最高的孕周
        return group.loc[group['达标比例'].idxmax(), '实际孕周']

    optimal_weeks = results.groupby('BMI分组').apply(find_optimal_week).reset_index()
    optimal_weeks.columns = ['BMI分组', '最佳孕周']

    # 5. 评估风险
    def assess_risk(week):
        if week <= 12:
            return '低风险'
        elif week <= 27:
            return '高风险'
        else:
            return '极高风险'

    optimal_weeks['风险等级'] = optimal_weeks['最佳孕周'].apply(assess_risk)

    # 6. 可视化结果
    plt.figure(figsize=(12, 6))
    for bmi_group in bmi_labels:
        group_data = results[results['BMI分组'] == bmi_group]
        plt.plot(group_data['实际孕周'], group_data['达标比例'], label=bmi_group, marker='o')

    plt.axhline(y=0.8, color='r', linestyle='--', label='80%达标阈值')
    plt.xlabel('孕周')
    plt.ylabel('Y染色体达标比例')
    plt.title('不同BMI分组在不同孕周的Y染色体达标比例')
    plt.legend()
    plt.grid()
    # plt.show()
    plt.savefig( os.path.join( out_dir , '不同BMI分组在不同孕周的Y染色体达标比例.png' ) )


    # 7. 输出最终建议
    print("\n=== 最佳BMI分组及检测时间建议 ===")
    print(optimal_weeks.sort_values('最佳孕周'))

    # 假设Y染色体浓度测量误差为±0.5%
    measurement_error = 0.005

    # 计算考虑误差后的达标比例
    male_df['达标下限'] = ((male_df['Y染色体浓度'] - measurement_error) >= 0.04).astype(int)
    male_df['达标上限'] = ((male_df['Y染色体浓度'] + measurement_error) >= 0.04).astype(int)

    # 重新计算分组结果（保守估计）
    results_error = male_df.groupby(['BMI分组', '实际孕周'])[['达标下限', '达标']].mean().reset_index()
    results_error.rename(columns={'达标下限': '保守达标比例', '达标': '原始达标比例'}, inplace=True)

    # 比较误差前后的差异
    print("\n=== 测量误差影响分析 ===")
    print("误差前 vs 误差后的达标比例变化:")
    print(results_error.groupby('BMI分组')[['原始达标比例', '保守达标比例']].mean())

# 考虑多因素 进行 分组
from sklearn.ensemble import RandomForestRegressor
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
def stage9_1():
    # 创建输出文件夹
    out_dir = "ques3_prepare_data2"
    os.makedirs(out_dir, exist_ok=True)

    df = pd.read_csv('ques1_prepare_data/sheet1_男胎检测数据_fillnan_sifted_outlier.csv')

    # 处理 df 中 怀孕次数 把 ≥3 全部换乘3
    # 将怀孕次数大于等于3的值全部替换为3
    df['怀孕次数'] = df['怀孕次数'].apply(lambda x: 3 if x == '≥3' else x)
    # df.to_csv('ques1_prepare_data/sheet1_男胎检测数据_fillnan_sifted_outlier.csv')

    # 数据预处理
    # 筛选男胎数据并计算实际孕周
    male_df = df[df['Y染色体浓度'] > 0].copy()
    male_df['末次月经'] = pd.to_datetime(male_df['末次月经'])
    male_df['检测日期'] = pd.to_datetime(male_df['检测日期'])
    male_df['实际孕周'] = (male_df['检测日期'] - male_df['末次月经']).dt.days / 7
    male_df['达标'] = (male_df['Y染色体浓度'] >= 0.04).astype(int)

    # 计算达标所需的最早孕周（每个孕妇）
    earliest_达标 = male_df[male_df['达标'] == 1].groupby('孕妇代码')['实际孕周'].min()
    male_df['达标孕周'] = male_df['孕妇代码'].map(earliest_达标)

    male_df.to_csv('ques1_prepare_data/sheet1_男胎检测数据_fillnan_sifted_outlier.csv')


    # 特征工程：考虑多种因素
    features = ['年龄', '身高', '体重', '孕妇BMI', '怀孕次数', '生产次数',
                '唯一比对的读段数  ','GC含量',
                '13号染色体的Z值','18号染色体的Z值','21号染色体的Z值',
                'X染色体的Z值',
                '13号染色体的GC含量','18号染色体的GC含量','21号染色体的GC含量']
    target = '达标孕周'

    # 相关性分析
    correlation_matrix = male_df[features + [target]].corr()
    plt.figure(figsize=(15, 15))
    sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm', center=0)
    plt.title('各因素与Y染色体达标孕周的相关系数矩阵')
    plt.savefig( os.path.join( out_dir , '各因素与Y染色体达标孕周的相关系数矩阵.png' ) )
    plt.close()
    # plt.show()

    # 特征重要性分析（使用随机森林）
    X = male_df[features].fillna(male_df[features].median())
    y = male_df[target].fillna(male_df[target].median())

    rf = RandomForestRegressor(n_estimators=100, random_state=42)
    rf.fit(X, y)

    # 可视化特征重要性
    feature_importance = pd.DataFrame({
        'feature': features,
        'importance': rf.feature_importances_
    }).sort_values('importance', ascending=False)

    plt.figure(figsize=(15, 6))
    plt.barh(feature_importance['feature'], feature_importance['importance'])
    plt.xlabel('特征重要性')
    plt.title('影响Y染色体达标时间的关键因素')
    plt.gca().invert_yaxis()
    # plt.show()

    plt.savefig( os.path.join( out_dir , '影响Y染色体达标时间的关键因素.png' ) )
    plt.close()

from sklearn.metrics import silhouette_score, calinski_harabasz_score
def stage9_2():
    # 创建输出文件夹
    out_dir = "ques3_prepare_data2"
    os.makedirs(out_dir, exist_ok=True)

    male_df = pd.read_csv('ques1_prepare_data/sheet1_男胎检测数据_fillnan_sifted_outlier.csv')

    # important_features = ['年龄', '身高', '体重', '孕妇BMI',
    #                         '唯一比对的读段数  ','GC含量',
    #                         '13号染色体的Z值','18号染色体的Z值','21号染色体的Z值',
    #                         'X染色体的Z值',
    #                         '13号染色体的GC含量','18号染色体的GC含量','21号染色体的GC含量']    
    important_features = ['体重',
                            '孕妇BMI',
                            '身高', 
                            '21号染色体的Z值',
                            '18号染色体的GC含量',
                            '年龄',  
                            'X染色体的Z值',
                            '13号染色体的Z值',
                            '唯一比对的读段数  ',
                            '21号染色体的GC含量','13号染色体的GC含量',
                            'GC含量', '18号染色体的Z值',
                            ]
    # 数据标准化
    scaler = StandardScaler()
    scaled_features = scaler.fit_transform(male_df[important_features].fillna(male_df[important_features].median()))

    # 使用肘部法则确定最佳聚类数量
    inertia = []
    silhouette_scores = []
    ch_scores = []
    center_range = range(2,20)
    for k in center_range:
        kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)
        kmeans.fit(scaled_features)
        inertia.append(kmeans.inertia_)
        # 修正：使用labels_而不是inertia_
        labels = kmeans.labels_
        silhouette_scores.append(silhouette_score(scaled_features, labels))
        ch_scores.append(calinski_harabasz_score(scaled_features, labels))
    


    # plt.figure(figsize=(10, 6))
    # plt.plot(center_range, inertia, marker='o')
    # plt.xlabel('聚类数量')
    # plt.ylabel(' inertia')
    # plt.title('肘部法则确定最佳聚类数量')
    # plt.show()
    # 绘制多个评估指标
    fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 5))

    # Inertia
    ax1.plot(center_range, inertia, marker='o')
    ax1.set_title('肘部法则')
    ax1.set_xlabel('聚类数量')
    ax1.set_ylabel('Inertia')

    # 轮廓系数（越接近1越好）
    ax2.plot(center_range, silhouette_scores, marker='o', color='green')
    ax2.set_title('轮廓系数')
    ax2.set_xlabel('聚类数量')
    ax2.set_ylabel('Silhouette Score')

    # Calinski-Harabasz指数（越高越好）
    ax3.plot(center_range, ch_scores, marker='o', color='red')
    ax3.set_title('CH指数')
    ax3.set_xlabel('聚类数量')
    ax3.set_ylabel('CH Score')

    plt.savefig( os.path.join( out_dir , '肘部法则确定最佳聚类数量.png' ) )
    plt.close()

    # 基于多因素聚类
    # optimal_clusters = 4  # 根据肘部法则确定
    optimal_clusters = silhouette_scores.index(max(silhouette_scores)) + 2
    print(f"自动选择的最佳聚类数: {optimal_clusters}")

    kmeans = KMeans(n_clusters=optimal_clusters, random_state=42)
    male_df['多因素聚类'] = kmeans.fit_predict(scaled_features)

    
    # 分析每个聚类的特征分布
    cluster_profiles = male_df.groupby('多因素聚类')[important_features + ['达标孕周']].agg(['mean', 'std'])
    print("各聚类群体的特征分布:")
    print(cluster_profiles)

    # 为每个多因素聚类群体分析最佳检测时间
    def analyze_optimal_time(cluster_data):
        # 计算每个孕周的达标比例
        week_range = np.arange(8, 25, 1)  # 8-24周
        results = []
        
        for week in week_range:
            # 计算在该孕周及之前的达标比例
            early_data = cluster_data[cluster_data['实际孕周'] <= week]
            if len(early_data) > 0:
                达标比例 = early_data['达标'].mean()
                results.append({'孕周': week, '达标比例': 达标比例, '样本量': len(early_data)})
        
        results_df = pd.DataFrame(results)
        
        # 找到达标比例≥80%的最早孕周
        optimal_week = results_df[results_df['达标比例'] >= 0.8]['孕周'].min()
        if pd.isna(optimal_week):
            optimal_week = results_df.loc[results_df['达标比例'].idxmax(), '孕周']
        
        return optimal_week, results_df

    # 为每个聚类群体确定最佳检测时间
    cluster_analysis = {}
    for cluster_id in range(optimal_clusters):
        cluster_data = male_df[male_df['多因素聚类'] == cluster_id]
        optimal_week, weekly_results = analyze_optimal_time(cluster_data)
        cluster_analysis[cluster_id] = {
            'optimal_week': optimal_week,
            'weekly_results': weekly_results,
            'cluster_size': len(cluster_data),
            'cluster_features': cluster_data[important_features].mean().to_dict()
        }

    # 可视化每个群体的达标曲线
    plt.figure(figsize=(12, 8))
    for cluster_id, analysis in cluster_analysis.items():
        plt.plot(analysis['weekly_results']['孕周'], 
                analysis['weekly_results']['达标比例'], 
                label=f'聚类{cluster_id}', 
                linewidth=2)
        plt.axvline(x=analysis['optimal_week'], linestyle='--', alpha=0.7)

    plt.axhline(y=0.8, color='red', linestyle=':', label='80%达标阈值')
    plt.axvline(x=12, color='orange', linestyle='--', label='早期风险边界')
    plt.xlabel('孕周')
    plt.ylabel('达标比例')
    plt.title('多因素聚类群体的Y染色体达标曲线')
    plt.legend()
    plt.grid(True)
    plt.savefig( os.path.join( out_dir , '多因素聚类群体的Y染色体达标曲线.png' ) )
    plt.close()


    # 定义风险函数
    def assess_risk(week):
        if week <= 12:
            return '低风险', 1
        elif week <= 20:
            return '中风险', 2
        elif week <= 27:
            return '高风险', 3
        else:
            return '极高风险', 4



    # 分析每个聚类的BMI范围
    cluster_bmi_ranges = male_df.groupby('多因素聚类')['孕妇BMI'].agg(['min', 'max', 'mean']).reset_index()
    cluster_bmi_ranges['BMI区间'] = cluster_bmi_ranges.apply(
        lambda x: f"{x['min']:.1f} - {x['max']:.1f}", axis=1
    )
    print("=== 各聚类的BMI区间分布 ===")
    print(cluster_bmi_ranges[['多因素聚类', 'min', 'max', 'mean', 'BMI区间']])

    # 重新定义群体特征描述
    final_recommendations = []
    for cluster_id in range(optimal_clusters):
        cluster_data = male_df[male_df['多因素聚类'] == cluster_id]
        
        # 获取该聚类的实际BMI范围
        bmi_min = cluster_data['孕妇BMI'].min()
        bmi_max = cluster_data['孕妇BMI'].max()
        bmi_mean = cluster_data['孕妇BMI'].mean()
        
        # 获取其他特征的范围
        age_min, age_max = cluster_data['年龄'].min(), cluster_data['年龄'].max()
        weight_min, weight_max = cluster_data['体重'].min(), cluster_data['体重'].max()
        
        optimal_week, weekly_results = analyze_optimal_time(cluster_data)
        risk_level, risk_score = assess_risk(optimal_week)
        
        recommendation = {
            '分组ID': cluster_id,
            'BMI区间': f"{bmi_min:.1f} - {bmi_max:.1f}",
            '平均BMI': f"{bmi_mean:.1f}",
            '年龄范围': f"{age_min:.0f} - {age_max:.0f}岁",
            '体重范围': f"{weight_min:.1f} - {weight_max:.1f}kg",
            '最佳检测孕周': optimal_week,
            '风险等级': risk_level,
            '样本数量': len(cluster_data),
            '12周达标比例': weekly_results[weekly_results['孕周'] <= 12]['达标比例'].max() if any(weekly_results['孕周'] <= 12) else 0
        }
        final_recommendations.append(recommendation)

    # 创建最终的分组建议表
    final_df = pd.DataFrame(final_recommendations)
    final_df['风险排序'] = final_df['最佳检测孕周'].rank()
    final_df = final_df.sort_values('风险排序')

    print("=== 基于多因素的最佳分组建议（含明确区间）===")
    print(final_df.to_string(index=False))

    # 可视化展示各组的BMI区间
    plt.figure(figsize=(12, 6))
    # colors = ['lightblue', 'lightgreen', 'lightcoral', 'gold']

    for i, (_, row) in enumerate(final_df.iterrows()):
        bmi_range = list(map(float, row['BMI区间'].split(' - ')))
        plt.barh(i, width=bmi_range[1]-bmi_range[0], left=bmi_range[0], 
                alpha=0.7, edgecolor='black')
        plt.text((bmi_range[0] + bmi_range[1])/2, i, 
                f"分组{i}\n{row['BMI区间']}", 
                ha='center', va='center', fontweight='bold')

    plt.xlabel('BMI值')
    plt.ylabel('分组')
    plt.title('各分组的BMI区间分布')
    plt.yticks(range(len(final_df)), [f'分组{i}' for i in range(len(final_df))])
    plt.grid(axis='x', linestyle='--', alpha=0.7)
    plt.tight_layout()
    # plt.show()
    plt.savefig( os.path.join( out_dir , '各分组的BMI区间分布.png' ) )
    plt.close()
    
    # 建议重新命名：
    group_names = {
        0: "高危肥胖组(BMI>30)",
        1: "中度超重组(BMI26-30)", 
        2: "高龄超重组(BMI27-36)",
        3: "标准超重组(BMI26-37)"
    }

    # 添加分组名称到最终建议
    final_df['分组名称'] = final_df['分组ID'].map(group_names)
    
    print("\n=== 临床分组建议 ===")
    for _, row in final_df.iterrows():
        print(f"\n{row['分组名称']}:")
        print(f"  • BMI范围: {row['BMI区间']}")
        print(f"  • 年龄分布: {row['年龄范围']}")
        print(f"  • 体重范围: {row['体重范围']}")
        print(f"  • 最佳检测时间: {row['最佳检测孕周']:.1f}周 ({row['风险等级']})")
        print(f"  • 12周预期达标率: {row['12周达标比例']:.1%}")
        print(f"  • 样本数量: {row['样本数量']}例")

    # 分析每个区间的边界特征
    print("\n=== 分组边界详细分析 ===")
    for i in range(optimal_clusters):
        cluster_data = male_df[male_df['多因素聚类'] == i]
        
        print(f"\n分组{i} ({group_names[i]}):")
        print(f"  BMI范围: {cluster_data['孕妇BMI'].min():.1f} - {cluster_data['孕妇BMI'].max():.1f}")
        print(f"  包含样本: {len(cluster_data)}例")
        print(f"  典型特征: {cluster_data['年龄'].mean():.1f}岁, {cluster_data['体重'].mean():.1f}kg")
        
        # 分析边界案例
        boundary_cases = cluster_data.nlargest(3, '孕妇BMI')
        print(f"  最高BMI案例: {boundary_cases['孕妇BMI'].values[0]:.1f} (孕周: {boundary_cases['实际孕周'].values[0]:.1f})")


def stage_9_2_optim():
    import os
    import numpy as np
    import pandas as pd
    import matplotlib.pyplot as plt

    from sklearn.model_selection import train_test_split
    from sklearn.metrics import roc_auc_score, accuracy_score
    from sklearn.ensemble import GradientBoostingClassifier
    from sklearn.cluster import KMeans
    from sklearn.metrics import silhouette_score
    from sklearn.preprocessing import StandardScaler
    from sklearn.tree import DecisionTreeRegressor, export_text

    # ========= 基本设置 =========
    out_dir = "ques3_prepare_data2"
    os.makedirs(out_dir, exist_ok=True)

    df = pd.read_csv('ques1_prepare_data/sheet1_男胎检测数据_fillnan_sifted_outlier.csv')
    df = df.copy()

    # -------- 列名鲁棒处理（有的表里“唯一比对的读段数”可能有空格） --------
    reads_col_candidates = ['唯一比对的读段数  ', '唯一比对的读段数', '唯一比对的读段数（唯一比对）']
    reads_col = next((c for c in reads_col_candidates if c in df.columns), None)

    # -------- 基础校验 --------
    required_cols = ['达标', '实际孕周', '年龄', '身高', '体重', '孕妇BMI',
                     'GC含量', '13号染色体的Z值', '18号染色体的Z值', '21号染色体的Z值', 'X染色体的Z值',
                     '13号染色体的GC含量', '18号染色体的GC含量', '21号染色体的GC含量']
    if reads_col:
        required_cols.append(reads_col)

    missing = [c for c in required_cols if c not in df.columns]
    if missing:
        raise ValueError(f"缺少必要列: {missing}")

    # 保证达标为0/1
    df['达标'] = (df['达标'] > 0).astype(int)

    # 读段数做log1p降尺度（若有）
    if reads_col:
        df[reads_col] = np.log1p(df[reads_col])

    # ========= 特征工程 =========
    base_feats = ['年龄', '身高', '体重', '孕妇BMI',
                  'GC含量', '13号染色体的Z值', '18号染色体的Z值', '21号染色体的Z值', 'X染色体的Z值',
                  '13号染色体的GC含量', '18号染色体的GC含量', '21号染色体的GC含量']
    if reads_col:
        base_feats.append(reads_col)

    # 与孕周的关键交互：周*BMI、周*年龄、周*读段数、周*21Z
    df['周_BMI'] = df['实际孕周'] * df['孕妇BMI']
    df['周_年龄'] = df['实际孕周'] * df['年龄']
    if reads_col:
        df['周_读段'] = df['实际孕周'] * df[reads_col]
    else:
        df['周_读段'] = 0.0
    df['周_Z21'] = df['实际孕周'] * df['21号染色体的Z值']

    model_feats = base_feats + ['实际孕周', '周_BMI', '周_年龄', '周_读段', '周_Z21']

    # 缺失值填充
    df[model_feats] = df[model_feats].fillna(df[model_feats].median())

    X = df[model_feats]
    y = df['达标']

    # ========= 建模：p(达标 | 周, 特征) =========
    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.2, random_state=42, stratify=y
    )

    clf = GradientBoostingClassifier(
        n_estimators=400,
        learning_rate=0.03,
        max_depth=3,
        subsample=0.9,
        random_state=42
    )
    clf.fit(X_train, y_train)

    # 简单验证
    y_proba = clf.predict_proba(X_test)[:, 1]
    auc = roc_auc_score(y_test, y_proba)
    acc = accuracy_score(y_test, (y_proba >= 0.5).astype(int))
    print(f"[模型] AUC={auc:.3f}, ACC={acc:.3f}")

    # ========= 生成个人“达标概率轨迹” =========
    week_grid = np.arange(8, 25, 1)  # 8~24周
    N = len(df)
    W = len(week_grid)
    traj = np.zeros((N, W), dtype=float)

    # 为了批量预测，构造每周的特征副本
    base_X_for_pred = df[base_feats].reset_index(drop=True)
    if reads_col is None:
        # 如果没有读段列，确保周_读段被创建为0
        pass

    for j, w in enumerate(week_grid):
        tmp = base_X_for_pred.copy()
        # 动态列
        tmp['实际孕周'] = w
        tmp['周_BMI'] = w * df['孕妇BMI'].values
        tmp['周_年龄'] = w * df['年龄'].values
        tmp['周_读段'] = w * (df[reads_col].values if reads_col else 0.0)
        tmp['周_Z21'] = w * df['21号染色体的Z值'].values

        tmp = tmp[model_feats].fillna(df[model_feats].median())
        traj[:, j] = clf.predict_proba(tmp)[:, 1]

    # ========= 由轨迹得到“个人最佳检测周”与关键指标 =========
    def earliest_week_reach(p, weeks, thr=0.8):
        idx = np.where(p >= thr)[0]
        if len(idx) == 0:
            # 若从未达阈值，取最大概率对应周
            return weeks[np.argmax(p)]
        return weeks[idx[0]]

    pred_opt_week = np.array([earliest_week_reach(traj[i], week_grid, 0.8) for i in range(N)])
    pred_rate12 = traj[:, np.where(week_grid <= 12)[0]].max(axis=1)  # 12周内最高达标概率
    slope_8_12 = traj[:, np.where(week_grid == 12)[0][0]] - traj[:, np.where(week_grid == 8)[0][0]]

    # ========= 按“曲线形状”聚类 =========
    # 标准化每个“周点”特征，避免某些周的方差支配
    scaler = StandardScaler()
    traj_std = scaler.fit_transform(traj)

    sil_scores = []
    k_candidates = list(range(2, 7))
    models = []
    for k in k_candidates:
        km = KMeans(n_clusters=k, random_state=42, n_init=10)
        labels = km.fit_predict(traj_std)
        score = silhouette_score(traj_std, labels)
        sil_scores.append(score)
        models.append((km, labels))

    best_idx = int(np.argmax(sil_scores))
    best_k = k_candidates[best_idx]

    # 若最优k=2，但k=3的分数非常接近（<=0.02差距），偏向3类以提升临床可用性
    if best_k == 2:
        idx_k2 = k_candidates.index(2)
        if 3 in k_candidates:
            idx_k3 = k_candidates.index(3)
            if sil_scores[idx_k2] - sil_scores[idx_k3] <= 0.02:
                best_k = 3
                best_idx = idx_k3

    km, labels = models[best_idx]
    print(f"[曲线聚类] 选择的簇数: {best_k}，Silhouette={sil_scores[best_idx]:.3f}")

    df['曲线聚类'] = labels
    df['预测最佳周'] = pred_opt_week
    df['预测12周内概率'] = pred_rate12
    df['8-12周概率提升'] = slope_8_12

    # ========= 组内统计与图 =========
    groups = []
    plt.figure(figsize=(10, 7))
    for g in range(best_k):
        idx = np.where(labels == g)[0]
        mean_curve = traj[idx].mean(axis=0)
        std_curve = traj[idx].std(axis=0)

        # 画组平均曲线
        plt.plot(week_grid, mean_curve, label=f'组{g} (n={len(idx)})', linewidth=2)
        # 阈值与参考线
        # （阴影可选，避免花哨：plt.fill_between(week_grid, np.maximum(0, mean_curve-std_curve), np.minimum(1, mean_curve+std_curve), alpha=0.15)）

        # 组汇总
        g_opt = pred_opt_week[idx].mean()
        g_rate12 = pred_rate12[idx].mean()
        g_slope = slope_8_12[idx].mean()
        groups.append({
            '组别': g,
            '样本数': int(len(idx)),
            '平均最佳检测周': round(float(g_opt), 2),
            '12周内平均达标概率': round(float(g_rate12), 3),
            '8-12周平均提升': round(float(g_slope), 3),
            'BMI均值': round(float(df.loc[idx, '孕妇BMI'].mean()), 2),
            '年龄均值': round(float(df.loc[idx, '年龄'].mean()), 2),
            '体重均值': round(float(df.loc[idx, '体重'].mean()), 2),
        })

    plt.axhline(0.8, linestyle='--', alpha=0.6)
    plt.xlabel('孕周'); plt.ylabel('达标概率'); plt.title('各曲线簇的平均达标概率曲线')
    plt.legend()
    plt.grid(alpha=0.3)
    plt.tight_layout()
    plt.savefig(os.path.join(out_dir, '曲线聚类_平均达标概率曲线.png'))
    plt.close()

    # 输出组表
    group_df = pd.DataFrame(groups).sort_values('平均最佳检测周')
    print("\n=== 按达标曲线分组的建议 ===")
    print(group_df.to_string(index=False))

    # ========= 风险分级（基于曲线而非静态特征） =========
    # 规则：
    # - 若组的 12周内平均达标概率 >=0.9 且 平均最佳周<=12 => 低风险（早检即可）
    # - 若 平均最佳周在(12,14] 或 12周概率在[0.75,0.9) => 中风险（建议12-14周）
    # - 其他 => 高风险（>14周）
    def risk_from_curve(avg_opt_week, avg_rate12):
        if avg_rate12 >= 0.9 and avg_opt_week <= 12:
            return '低风险'
        elif (12 < avg_opt_week <= 14) or (0.75 <= avg_rate12 < 0.9):
            return '中风险'
        else:
            return '高风险'

    group_df['风险等级'] = group_df.apply(
        lambda r: risk_from_curve(r['平均最佳检测周'], r['12周内平均达标概率']), axis=1
    )
    print("\n=== 分组与风险等级 ===")
    print(group_df[['组别','样本数','平均最佳检测周','12周内平均达标概率','风险等级']].to_string(index=False))

    # 保存表
    group_df.to_csv(os.path.join(out_dir, '曲线聚类_分组建议.csv'), index=False)

    # ========= 可解释性：用小决策树解释“预测最佳周” =========
    # 用原始可控特征（不含“实际孕周”），回归“预测最佳周”，导出规则
    explain_feats = base_feats.copy()
    if '孕妇BMI' not in explain_feats:
        explain_feats.append('孕妇BMI')

    X_explain = df[explain_feats].fillna(df[explain_feats].median())
    y_explain = df['预测最佳周'].values

    tree = DecisionTreeRegressor(max_depth=3, min_samples_leaf=50, random_state=42)
    tree.fit(X_explain, y_explain)

    rules = export_text(tree, feature_names=list(X_explain.columns), decimals=2)
    with open(os.path.join(out_dir, '解释_最佳检测周_规则.txt'), 'w', encoding='utf-8') as f:
        f.write(rules)

    print("\n=== 可解释规则（预测最佳检测周） ===")
    print(rules)

    # ========= 最终一句话总结 =========
    print("\n【结论】已基于“达标概率-孕周曲线”完成分组。请查看：")
    print(f" - {os.path.join(out_dir, '曲线聚类_平均达标概率曲线.png')}")
    print(f" - {os.path.join(out_dir, '曲线聚类_分组建议.csv')}")
    print(f" - {os.path.join(out_dir, '解释_最佳检测周_规则.txt')}")


if __name__=='__main__':
    main()